home *** CD-ROM | disk | FTP | other *** search
- /******************************************************************************
- * BspCoxDB.c - Bspline evaluation using Cox - de Boor recursive algorithm. *
- *******************************************************************************
- * Written by Gershon Elber, Aug. 90. *
- ******************************************************************************/
-
- #include <ctype.h>
- #include <stdio.h>
- #include <string.h>
- #include "cagd_loc.h"
-
- /*****************************************************************************
- * DESCRIPTION: M
- * Returns a pointer to a static data, holding the value of the curve at M
- * the prescribed parametric location t. M
- * Uses the recursive Cox de-Boor algorithm, to evaluate the spline, which M
- * is not very efficient if many evaluations of the same curve are necessary M
- * Use knot insertion when multiple evaluations are to be performed. M
- * *
- * PARAMETERS: M
- * Crv: To evaluate at the given parametric location t. M
- * t: The parameter value at which the curve Crv is to be evaluated. M
- * *
- * RETURN VALUE: M
- * CagdRType *: A vector holding all the coefficients of all components M
- * of curve Crv's point type. If for example the curve's M
- * point type is P2, the W, X, and Y will be saved in the M
- * first three locations of the returned vector. The first M
- * location (index 0) of the returned vector is reserved for M
- * the rational coefficient W and XYZ always starts at second M
- * location of the returned vector (index 1). M
- * *
- * KEYWORDS: M
- * BspCrvEvalCoxDeBoor, evaluation, Bsplines M
- *****************************************************************************/
- CagdRType *BspCrvEvalCoxDeBoor(CagdCrvStruct *Crv, CagdRType t)
- {
- static CagdRType Pt[CAGD_MAX_PT_COORD];
- CagdBType
- IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
- CagdRType *p, *BasisFunc;
- int i, j, l, IndexFirst,
- k = Crv -> Order,
- Length = Crv -> Length,
- MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
-
- BasisFunc = BspCrvCoxDeBoorBasis(Crv -> KnotVector, k,
- CAGD_CRV_PT_LST_LEN(Crv),
- t, &IndexFirst);
-
- /* And finally multiply the basis functions with the control polygon. */
- for (i = IsNotRational; i <= MaxCoord; i++) {
- Pt[i] = 0;
- p = Crv -> Points[i];
- for (j = IndexFirst, l = 0; l < k; )
- Pt[i] += p[j++ % Length] * BasisFunc[l++];
- }
-
- return Pt;
- }
-
- /******************************************************************************
- * DESCRIPTION: M
- * Returns a pointer to a vector of size Order, holding values of the non M
- * zero basis functions of a given curve at given parametric location t. M
- * This vector SHOULD NOT BE FREED. Although it is dynamically allocated, M
- * the returned pointer does not point to the beginning of this memory and it M
- * it be maintained by this routine (i.e. it might be freed next time this M
- * routine is being called). M
- * IndexFirst returns the index of first non zero basis function for the M
- * given parameter value t. M
- * Uses the recursive Cox de Boor algorithm, to evaluate the Bspline basis M
- * functions. M
- * Algorithm: M
- * Use the following recursion relation with B(i,0) == 1. M
- * M
- * t - t(i) t(i+k) - t V
- * B(i,k) = --------------- B(i,k-1) + --------------- B(i+1,k-1) V
- * t(i+k-1) - t(i) t(i+k) - t(i+1) V
- * M
- * Starting with constant Bspline (k == 0) only one basis function is non M
- * zero and is equal to one. This is the constant Bspline spanning interval M
- * t(i)...t(i+1) such that t(i) <= t < t(i+1). We then raise this constant M
- * Bspline to the prescribed Order and find in this process all the basis M
- * functions that are non zero in t for order Order. Sound simple hah!? M
- * *
- * PARAMETERS: M
- * KnotVector: To evaluate the Bspline Basis functions for this space. M
- * Order: Of the geometry. M
- * Len: Number of control points in the geometry. The length of M
- * KnotVector is equal to Len + Order. M
- * t: At which the Bspline basis functions are to be evaluated. M
- * IndexFirst: Index of the first Bspline basis function that might be M
- * non zero. M
- * *
- * RETURN VALUE: M
- * CagdRType *: A vector of length Order thats holds the values of the M
- * Bspline basis functions for the given t. A Bspline of M
- * order Order might have at most Order non zero basis M
- * functions that will hence start at IndexFirst and upto M
- * (*IndexFirst + Order - 1). M
- * *
- * KEYWORDS: M
- * BspCrvCoxDeBoorBasis, evaluation, Bsplines M
- *****************************************************************************/
- CagdRType *BspCrvCoxDeBoorBasis(CagdRType *KnotVector,
- int Order,
- int Len,
- CagdRType t,
- int *IndexFirst)
- {
- static CagdRType
- *B = NULL;
- CagdRType s1, s2, *BasisFunc;
- int i, l, Index,
- KVLen = Order + Len;
-
- if (!BspKnotParamInDomain(KnotVector, Len, Order, FALSE, t))
- CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
- if (t == KnotVector[Len])
- t -= IRIT_EPSILON;
- Index = BspKnotLastIndexLE(KnotVector, KVLen, t);
-
- /* Starting the recursion from constant splines - one spline is non */
- /* zero and is equal to one. This is the spline that starts at Index. */
- /* As We are going to reference index -1 we increment the buffer by one */
- /* and save 0.0 at index -1. We then initialize the constant spline */
- /* values - all are zero but the one from t(i) to t(i+1). */
- if (B != NULL)
- IritFree((VoidPtr) B);
- BasisFunc = B = (CagdRType *) IritMalloc(sizeof(CagdRType) * (1 + Order));
- *BasisFunc++ = 0.0;
- if (Index >= Len + Order - 1) {
- /* We are at the end of the parametric domain and this is open */
- /* end condition - simply return last point. */
- for (i = 0; i < Order; i++)
- BasisFunc[i] = (CagdRType) (i == Order - 1);
-
- *IndexFirst = Len - Order;
- return BasisFunc;
- }
- else
- for (i = 0; i < Order; i++)
- BasisFunc[i] = (CagdRType) (i == 0);
-
- /* Here is the tricky part. we raise these constant splines to the */
- /* required order of the curve Crv for the given parameter t. There are */
- /* at most order non zero function at param. value t. These functions */
- /* start at Index-order+1 up to Index (order functions overwhole). */
- for (i = 2; i <= Order; i++) { /* Goes through all orders... */
- for (l = i - 1; l >= 0; l--) { /* And all basis funcs. of order i. */
- s1 = (KnotVector[Index + l] - KnotVector[Index + l - i + 1]);
- s1 = APX_EQ(s1, 0.0) ? 0.0 : (t - KnotVector[Index + l - i + 1]) / s1;
- s2 = (KnotVector[Index + l + 1] - KnotVector[Index + l - i + 2]);
- s2 = APX_EQ(s2, 0.0) ? 0.0 : (KnotVector[Index + l + 1] - t) / s2;
-
- BasisFunc[l] = s1 * BasisFunc[l - 1] + s2 * BasisFunc[l];
- }
- }
-
- *IndexFirst = Index - Order + 1;
- return BasisFunc;
- }
-